Topological Entanglement Entropy in the Quantum Dimer Model 

on the Triangular Lattice 
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A characterization of topological order in terms of bi-partite entanglement was proposed recently 
[A. Kitaev and J. Preskill, Phys. Rev. Lett. 96, 110404 (2006); M. Levin and X.-G. Wen, ibid, 
110405]. It was argued that in a topological phase there is a universal additive constant in the 
entanglement entropy, called the topological entanglement entropy, which reflects the underlying 
gauge theory for the topological order. In the present paper, we evaluate numerically the topological 
entanglement entropy in the ground-states of a quantum dimer model on the triangular lattice, 
which is known to have a dimer liquid phase with Z2 topological order. We examine the two original 
constructions to measure the topological entropy by combining entropies on plural areas, and we 
observe that in the large-area limit they both approach the value expected for Z2 topological order. 
We also consider the entanglement entropy on a topologically non-trivial "zigzag" area and propose 
to use it as another way to measure the topological entropy. 

PACS numbers: 75.10.Jm, 03.65.Ud, 05.30.-d 



I. INTRODUCTION 

Exotic phenomena in quantum many-body systems are 
accompanied by non-trivial patterns of entanglement in 
ground-state wave-functions. One useful measure of en- 
tanglement for a many-body state jvP) is the entangle- 
ment entropy Sq between a part £1 of the system and the 
rest of the system, ft. It is defined as the von Neumann 
entropy of the reduced density matrix pq obtained by 
tracing out the degrees of freedom of f2: 



Sn = -Trpnlnpn, 



pn = Tr|*)(*| 

Q 



(1) 



It has been clarified in the past few years that some 
important properties of a quantum ground-state are en- 
coded in the size-dependence of Sq. For a system with 
short-range correlations only, £1 and f2 correlate only in 
the vicinity of the boundary separating them and thus 
the entanglement entropy, scales with the size of the 
boundary (boundary law) .El However, at a critical point 
with algebraically decaying correlations, the scaling of en- 
tanglement entropy exhibits a universal logarithmic cor- 
rection characterizing the criticality. Specifically, in a 
one-dimensional quantum critical system described by a 
conformal field theory (CFT), the entanglement entropy 
shows a logarithmic scaling law with a coefficient deter- 
mined by the central charge of the CFT.cl In some two- 
dimensional quantum critical states, the entanglement 
entropy also contains a universal contribution, related to 
the geometry of the subsystem.El 

Another type of non-trivial entanglement can exist in a 
system with topological orderuu Such a system exhibits 
degenerate ground-states separated from excited states 
by an energy gap, and this degeneracy, which depends on 
the topology of the entire system, cannot be ascribed to 
any type of conventional spontaneous symmetry break- 



ing. Indeed, it has been demonstrated in some models 
that these degenerate ground states cannot, be distin- 
guished by any local observable EBu PreskilB suggested 
that this degeneracy can be regarded as a global encoding 
of information reminiscent of quantum error-correcting 
codes and is a consequence of some long-distance entan- 
glement. A characterization of this global entanglement 
was realized recently by Kitaev and Preskill (KP)EJ and 
by Levin and Wen (LW)EJ. It was argued that, if O is a 
disk (in a two-dimensional system) with a smooth bound- 
ary of length L, the entanglement entropy scales as 



Sn = aL — 7 + 



(2) 



where the ellipsis represents terms which are negligible 
in the limit L — > 00. If the area f2 is not a disk and 
has m disconnected boundaries, the topological term —7 
in Eqn. (||) is multiplied by m. While the coefficient a 
depends on the microscopic details of the system, 7 is 
a universal constant characterizing topological order and 
was dubbed the topological entanglement entropy. In- 
deed, 7 measures the so-called total quantum dimension 
V of topological order by 7 = In T>. In the case of topo- 
logical order described by a discrete Abelian gauge theory 
(e.g., Z ra ), T> is equal to the number of elements in the 
gauge group. In general, it is difficult to separate the 
topological term —7 from the boundary term in Eqn. (Q) 
because, on a lattice, the discrete nature of the bound- 
ary makes it difficult to define unambiguously the length 
L. To solve this, KP and LW found some ways to define 
7 by forming a linear combination of the entanglement 
entropies on plural areas sharing some boundaries, and 
cancelling the boundary terms out to leave the topologi- 
cal term. KP and LW illustrated this idea using effective 
field theories and exactly solvable models. 

In this paper, we analyze the entanglement entropy 
in the quantum dimer model (QDM) on the triangular 
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lattic and examine the effectiveness of the proposal in 
numerical calculations of finite-size systems. This model 
is known to exhibit a dimer liquid phase with Z2 topp=. 
logical order in a finite interval in the parameter space. Li 
We mainly consider the Rokhsar-Kivelson (RK) point ,E2 
where the ground-states are exactly known and where 
the calculation of reduced density matrices (and thus en- 
tanglement entropy) amounts to counting the number of 
dimer coverings of the lattice satisfying some particu- 
lar constraints. We calculate the topological entangle- 
ment entropy numerically, and compare the result with 
7 = In 2 expected for Z 2 topological order. 

We .-comment on related systems here. Kitaev's 
modeE.3 is known to be the simplest solvable model with 
Z2 topological order, and the entanglement entropy of 
this model has been analyzed rigorously in Refs. [ll| , p^| 
and the value 7 = In 2 for the topological entropy was 
confirmed. The solvable QDM (kagome lattice) of Ref. II 
can be mapped onto Kitaev's model on the honeycomb 
lattice, and thus its entanglement entropy can be ana- 
lyzed in the same way. These models give elegant re- 
sults, but are too ideal for discussing generic features of 
topological order because they have a strictly zero spin- 
spin (or dimer-dimer in the QDM) correlation length and 
are completely free of finite-size effects. In this sense, our 
analysis on the QDM on the triangular lattice is a step to- 
ward more realistic systems - though we mainly consider 
the exact RK ground-states, they have a finite dimer- 
dimer correlation length and finite-size effects arise. In 
the same spirit but for another kind of topological order, 
the entanglement entropy of Laughlin wave functions was 



analyzed numerically in Ref. 17 

The paper is organized as follows. In Sec. ||, we give 
the basic definitions and settings in our analysis. In 



Sec. [II, we numerically analyze the properties of entan- 
glement entropy in the QDM on the triangular lattice. 
Especially, we examine the two constructions of topo- 
logical entanglement entropy proposed by KP and LW. 
Furthermore, we consider the entanglement entropy on 
a particular topologically non-trivial area and design an- 
other procedure to extract 7, which, for QDMs, turns out 
to give an accurate value even in relatively small systems. 
We then conclude in Sec. 



IV 



II. DEFINITIONS AND SETTINGS 



Model 



where the sum runs over all rhombi consisting of two 
neighbouring triangles and we set t > 0. At the Rokhsar- 
Kivelson (RK) point v = t, a ground-state is given ex- 
actly by the equal-amplitude superposition of all the 
dimer coverings £3 



|RJK) = 



1 



Ei c >> 



(4) 



where £ denotes the set of all the dimer coverings. This 
wave function e xhibits , exponentially-decaying dimer- 
dimer correlations! 1 tHM and is an example of liquid with 
no broken symmetries. 

This wave function is not the unique ground state if 
the lattice has a non-trivial topology (cylinder, torus, 
etc.). Let us focus on the case of the torus hereafter. 
We draw two incontractible loops Ai and A2 which pass 
through the bonds and wind around the torus in x and 
y directions respectively as in Fig. [p. We classify £ into 

four sets £ v with p = ++,+ — ,— + , , depending on 

the parity of the number of dimers crossing Ai and A2. 
The resultant sets £ p , called topological sectors, are not 
mixed by any local dimer move (and thus by any term in 
the Hamiltonian) . The spectrum of the Hamiltonian can 
therefore be determined separately in each sector. At the 
RK point, the ground-state in each sector is given by 



|RK;p) = 




(5) 



All these states have zero energy for the Hamiltonian (j^) 
and span a four-dimensional ground-state manifold. It 
has been shown analytically and numerically that the de- 
generacy of the ground states and the exponential decay 
of the dimer-dimer correlation at the RK point persist 
in a finite range in the parameter space, forming a liq- 
uid_|phasc-.with gapped excitations in 0.82(3) < v/t < 
Decreasing v/t further, the model enters a 
valence bond crystal (VBC) phasejffiith a large unit cell 
(12 sites), called VT2 x Vl2 VBC. 00 

The ground-state degeneracy in the liquid phase indi- 
cates that this phase is topologically ordered. It is indeed 
a realization of the deconfined phase of a Z2 (Ising) gauge 
theory,E3 where the requirement that physical states must 
be invariant under gauge transformations is played by 
the dimer hard-core constraint and where the role of the 
gauge flux piercing a plaquette is.pl a yed by a dimer-move 
operator around this plaquetteES'E^I The four ground- 
states correspond to the four possible choices to put (or 
not to put) a vortex through the two holes of the torus. 



We consider the QOMi on the triangular lattice defined 
by the HamiltoniamMliia 



*= E 



rhombi 



h.c. 



// // 



(3) 



B. Lattice 

The lattice is put on a torus and is defined by two 
vectors Ti and T2 specifying the periodicity. We mostly 
use lattices which are symmetric under 120° rotation, by 
setting 



Ti = lu + mv, 



(I + m)v, 



(6) 
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where I and to are integers and u and v are unit vec- 
tors as shown in Fig. |]. The total number of sites is 
given by N — I 2 + Im + to 2 . The lattices we consider 
have N = 16, 28, 36, 48, 52, 64, which correspond respec- 
tively tojX to) = (4, 0) , (4, 2) , (6, 0) , (4, 4) , (6, 2) , (8, 0) . 
In Sec. [II C , N = 100 (corresponding to (10,0)) is also 
studied. 




FIG. 1: Triangular lattices with periodic boundary condi- 
tions. 



C. Reduced density matrix 

To define a reduced density matrix (RDM) for the 
QDM, we must specify the local degrees of freedom of 
the model. To this end, we assign an Ising variable 
to each bond k of the lattice as in Ref. |2^ and identify 
the presence/ absence of a dimer on the bond as <7fc = +1 
and —1, respectively. Any physical configuration {<Jk} 
must satisfy the hard-core constraints: for each site of 
the lattice, there must be exactly one bond with a\~ = 1 
emanating from it. An area fl is defined as a set of bonds. 
We define the matrix element of the RDM of a ground- 
state I'F) as 



(ci\pn\c 2 ) = ^(c 1 ,c|1')(1'|c2,c), 



(7) 



where ci and C2 are dimer configurations on £1 and the 
sum is over all the dimer configurations c on O. Note that 
we set (c, c\^) =0 if (c, c) is an unphysical configuration 
(violating the hard-core constraint). 

Since the liquid phase under consideration exhibits de- 
generate ground states, we must specify for which state in 
the ground-state manifold we calculate the entanglement 
entropy. However, as long as the area is local, it was nu- 
merically demonstrated that the RDMs are identical for 
all states in the ground-state manifold, up to a correction 
which decays exponentially with the system sizc.0 Thus 
in this case we can take any state in the ground-state 
manifold. At the RK point, which we mainly consider 
in the following, we simply take the "equal-amplitude" 
state (Q) . The RDM of the "equal-amplitude" state can 



be calculated in a way described in the Appendix [A|, ei- 
ther by direct enumeration, or using Pfaffians. 



III. NUMERICAL RESULTS 

Here we present our numerical results. The idea of 
KP and LW should apply to the dimer liquid phase in 
0.82(3) < v/t < 1. The topological entropy for this phase 
is expected to be 7 = In 2 ~ 0.6931, reflecting Z2 topo- 
logical order. We mainly consider the RK point v/t = 1 
with exact ground-states (||) or (||) , and calculate the en- 
tanglement entropies using the methods in Appendix |a|. 
For v/t < 1, we perform Lanczos diagonalization of the 
Hamiltonian (||) for small systems (up to N = 36), and 
calculate the entanglement entropies in the ground-state. 



A. Circular areas 

We first consider the entanglement entropy on disks 
(areas with no holes) and discuss how the entanglement 
entropy scales with the extension of the area. Calcula- 
tions were done for the RK wave function (|[). As the 
choice of the area f2, we define circular areas in the fol- 
lowing way: we draw a circle with a radius R centred 
at a site or at an interior of a triangle and regard every 
bond whose midpoint is in the circle as an element of the 
area; see Fig. pi This definition causes an unavoidable 
ambiguity in the radius R — different radii can result 
in the same area. For example, the possible radius for 
the smallest site-centered area (consisting of six bonds) 

ranges in i? m i n = 0.5<i?<-^ = i? max . Here we ana- 
lyze the data taking this ambiguity into account. 

In Fig. ||, the values of Sq on circular areas are plotted 
versus the radius R. The different symbols correspond to 
different system sizes (from N = 16 to 52) and an hor- 
izontal bar specifies the interval [i? m i n , if m ax] (N = 52 
data points only, for clarity) The data from different sys- 
tem sizes almost coincide, showing the smallness of the 
finite-size effects. We fit the data for N = 52 by a linear 
relation using i? m i n or i? max We observe a rough agree- 
ment with the linear fitting in both cases as expected 
from the scaling form (||) . The lines intersect the vertical 
axis around —0.1 and —1.8 when using i? m i n and i? max , 
respectively. These values sandwich the expected value 
— In 2 ~ —0.6931 but are both away from it. We also 
fitted the data separately for site-centered and triangle- 
centered cases (not shown in the figure), but no essential 
difference was observed. These results show that a direct 
check of the scaling @ is difficult. 

In general, on a lattice, the boundary of ft is made of 
segments. If the sum of the segments is long enough, they 
contribute to the entanglement entropy by an amount 
proportional to the length. But in addition, we have 
to take into account the contribution coming from local 
correlations (between the regions fl and Cl) taking place 
in the vicinity of the angles between successive segments. 



If n is large, the contribution from these angles may be 
small (of order O(L ), compared to the boundary length 
L), but this contribution will still be of the same order 
as the topological term we are looking for. 

In the present case (circular areas), this ambiguity in 
defining a boundary length on the lattice appears as an 
ambiguity in the definition of R. To compute 7 in a well- 
defined way, we need to turn to the constructions using 
plural areas, which we discuss in the next subsection. 



(a) S00 



(b) S30< 




FIG. 2: Circular areas centered at (a) a site or (b) an interior 
of a triangle. As examples, areas with R = 2.5 and R — 2.47 
are shaded for (a) and (b), respectively. 




FIG. 4: Divisions of circular areas for the Kitaev-Preskill con- 
struction, (a) and (b): site-centered, R = 2.78. (c) and (d): 
triangle-centered, R = 2.84. 




_2 t- 1 1 1 1 1 L_ 

0.5 1 1.5 2 2.5 3 
Radius R 



FIG. 3: (color online) Entanglement entropy on circular ar- 
eas with radii R at the RK point. The ambiguity in R is 
indicated by horizontal bars (only for N = 52). The data for 
TV = 52 are fitted by lines using minimum or maximum radii. 
The resultant linear functions shown in the figure contains 
some numbers enclosed in parentheses, which indicates the 
standard errors in the last displayed digits. 



B. Construction of the topological entropy using 
plural areas 

KP and LW proposed two ways to extract the topo- 
logical constant Tijpdependently of the definition of the 
boundary length.t2llill The idea is to evaluate 7 by form- 
ing an appropriate linear combination of the entangle- 



ment entropies of different areas, so that the boundary 
contributions cancel out. 



1. Kitaev-Preskill construction 

In the KP construction^ we consider a circle and di- 
vide it into three "fans", A, B, and C. Then we form a 
linear combination 

■Stopo = Sa + Sb + Sc — Sab — Sbc — &ca + Sabc, (8) 

where Sxy— denotes the entanglement entropy on a com- 
posite area X U Y U • • • . In this combination, all the 
boundary contributions cancel out and a topological term 
—7 should remain. For example, let us consider the line 
separating A and B. The boundary contributions along 
this line appears in Sa and Sb with a plus sign and 
in —Sbc an d — Sca with a minus sign. Some attention 
should be paid to the triple point, in the vicinity of which 
the areas have different shapes and thus possibly different 
local contributions. Three areas form a 120 degree angle: 
A, B and C; three areas form a 240-degrees angle: BC, 
AC and AB. However, recalling that the entanglement 
entropies of an area and its complement are the same, the 
entropy of BC is equal to that of the complement of BC, 
which has the same shape with A in the vicinity of the 
triple point. Thus the local contributions from A and BC 
in the vicinity of the triple point should match. The same 
argument applies to every line and every corner, giving 
a cancellation of all the boundary contributions in Eqn. 
(||). Assuming the scaling (||), we expect S^ Q = —7 (for 
a large enough radius). 

We apply this idea to the present model. We di- 
vide a circle by three lines emanating from the center 
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to 



-0.2 



-0.4 



-0.6 



-0.8 



Kitaev-Preskill, RK point 

N=36 v 

N=48 A 

N=52 □ 

N=64 o 




R 



FIG. 5: (color online) Topological entanglement entropy from 
the Kitaev-Preskill construction (H) at the RK point. Exam- 
ples of areas are shown in Fig. Some explicit values for 
large radii are shown in Table | 



TABLE I: Some values of — Sjopo at the RK point, divided by 
the expected value In 2. The data for large radii are shown, 
and excellent agreement with the expectation can be seen. 



S00 case 
Radius R -£££ /ln2 



T30 case 
Radius R -S^ /ln2 





N = 52 


N = 64 




N = 52 


N = 64 


2.18 


0.9143 


0.9143 


2.57 


0.9291 


0.9283 


2.29 


0.9839 


0.9835 


2.75 


0.9618 


0.9513 


2.50 


0.9822 


0.9822 


2.84 


0.9965 


0.9518 


2.60 


0.9765 


0.9760 


2.93 


1.0910 


0.9635 


2.78 


1.0014 


0.9897 


3.01 


1.0910 


0.9635 


3.04 


1.3252 


0.9967 


3.18 




0.9649 


3.12 




0.9967 


3.25 




0.9898 




v/t 



FIG. 6: (color online) Kitaev-Preskill topological entropy (J]) 
as a function of v/t for TV = 36. In the large- J? limit, 
expected jump from In 2 in Z2 liquid phase 0.82(3) < v/t < 1, 
to some positive value in the VBC phase v/t < 0.82(3). 



as in Fig. |^. These lines are placed at angles 60 — 0, 
9 a + 120° - and 6 + 240° - measured from the (refer- 
ence) u direction. Here "—0" represents an infinitesimal 
shift for avoiding collisions between the points (midpoints 
of bonds) and the boundaries. For example, points at 
an angle #0 belong to A, not to C. We take 9o = 0° 
or 30° for site-centered circles (referred to as "S00" and 
"S30") and 9 = 30° or 90° for triangle-centered circles 
("T30" and "T90"). In these settings, the parts A,B,C 
are equivalent under 120° rotation, and we thus only need 
to calculate S^po = 3Sa — 3Sab + Sabc- 

We first consider the case of the RK wave function 
(U). In Fig. H the data of are plotted versus the 

radii R of the circles. As in the case of circular areas pre- 
sented in Fig. |[ finite-size effects are very small - except 
for the case where the circle ABC occupies a substantial 
part of the system, the data from different AT's almost 
coincide. In the largest system A" = 64, we can regard 
the data up to R < 3.1 as good approximation to the 
values in the infinite system. In all the cases, S^po de- 
creases almost monotonically with R and for large radii 
(specifically, 2.2 < R < 3.1) shows values which are very 



close to — In 2 (see Table |) , the expected value for a Z2 
topologically ordered state. 

Next we consider the region v/t < 1 of the Hamilto- 
nian (§). In Z 2 liquid phase 0.82(3) < v/t < 1, S££ Q 
is expected to show — In 2 in the large- R limit. On the 
other hand, in \/l2 x ^12 VBC phase v/t < 0.82(3), 
where discrete symmetries are spontaneously broken, the 
finite-size ground-state can be approximated by a linear 
superposition of 12-fold symmetry-broken states. In such 
a state, we conjecture that the entanglement entropy on 
a disk scales as Sq ~ aL + In d in the large- area limit, 
where d is the ground-state degeneracy and is equal to 12 
in the present case. The constant term hid is not topo- 
logical in the sense that the same value would appear 
even if fl had another geometry, unlike —7 in Eqn. (||). 
Note also that this constant is positive, in contrast to the 
negative topological term —7. Assuming this, the combi- 
nation (pi) should give In d in a symmetry-broken phase. 
Thus, S'Jopo i s expected to jump from a negative (topo- 
logical) value — In 2 to a positive (non-topological) value 
In 12 along with the transition from the liquid phase to 
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the VBC phase. Weperformed Lanczos diagonalization 
of the Hamiltonian (Q) for a lattice with N = 36 (which 
is the maximum size in our exact diagonalization calcu- 
lation and is compatible with \/l2 x Vl2 VBC ordering), 
and calculated in the ground-state, which lies in 

the sector p = in both the VBC and liquid phases 

on this lattice. The results are shown for two types of 
areas ("S00" and "T30") in Fig. |. Because the system 
and area sizes are rather small, we do not observe a jump 
at the transition. However, we can already observe some 
tendency: for fixed v/t, tends to decrease as a func- 
tion of R in the liquid side while it tends to increase in 
the VBC side. Some positive values of S^ Q in the VBC 
phase are also seen in "T30" case. 

2. Levin-Wen construction 

In the LW construction^, we consider an annulus di- 
vided into four pieces as in Fig. ^, and form a combina- 
tion 

= Sabcd — Sabc — Scda + Sac- (9) 

This combination is guaranteed to be non-positive from 
the strong, subadditivity inequality of entanglement 
entropies,Cj namely, 

&top° ~ Sxuy — Sx — Sy + SxnY < 0, (10) 

where X = AUBUCand Y = CUDUi. The combina- 
tion (ph is expected to give —27 for a topological phase 
and zero for a conventional phase (disordered, or with 
some symmetry-breaking order). 

In Fig. 0, an annulus is divided by four lines at angles 
6»o - 0, O + 60° + 0, O + 180° -O,0 O + 240° + 0. We 
consider only site-centered annuli, and we set 9q = 0° 
or 30° (again referred to as "S00" and "S30"). The re- 
sult for the RK wave function is shown in Fig. |^. R ln 
and i? ou t denote the inner and outer radii of the annu- 
lus respectively, and S™ Q 's are plotted as a function 
of .Rout- Up to i? out < 3.1, where the data for N = 64 
well approximate the values in the infinite system, we 
observe that S™ Q monotonically decreases with R out 
and approaches —2 In 2. Unfortunately, the convergence 
to —2 In 2 is not very clear up to this radius. In the 
LW construction, the requirement for the convergence 
is £ << Rin,Rout — Rin,L — 2i? out , where £ is the cor- 
relation length (~ 1 at RK point) and L — VN is the 
linear system size (or equivalently, the maximum possible 
2-Rout)- Thus, the LW construction suffers from stronger 
finite-area (not finite- N) effects than the KP construction 
which just requires £ << R, L — 2R. 

C. Zigzag area 

We design a different way to evaluate 7 using a thin 
"zigzag" area il winding around the torus as in Fig. 0. 




FIG. 7: Division of annular areas for the Levin- Wen construc- 
tion. 



5 -1 

_J 

O 
Q. 

c/f -1.2 



I I I I I 
(a) Levin-Wen, S00, RK point 


1 1 

N=36 v 
N=48 A 
N=52 □ 
N=64 O 


-2 In 2=- 1 .386 \\ 

i i i i i 


a \vi-32 

q R in =0.87 - 
1 1 



2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 
R out 



-0.6 




i i i i i 
I (b) Levin-Wen, S30, RK point 


1 1 

N=36 v 
N=48 A 
N=52 □ 


-0.8 






N=64 O 
R in =1.50 


-1 








-1.2 
-1.4 






\\)R irl =0.87 




-2 In 2=- 1 .386 \\ 


®R in =1.32 


-1.6 




1 1 1 1 1 


1 T 



2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 
R out 



FIG. 8: Topological entanglement entropy from the Levin- 
Wen construction. 



This area is invariant by translation in the x direction 
and all points (black circles in Fig. ^) are equivalent by 
symmetry. In contrast to the more complicated areas 
considered before, we expect the boundary (i.e., non- 
topological) contribution to Sq to be precisely propor- 
tional to l x , when l x is sufficiently larger than the correla- 
tion length £. In this new geometry, the thermodynamic 
behavior is obtained as soon as £ <C Ixi^yi whereas the 
KP construction requires £ << R,L — 2R, which is dif- 
ficult to reach in exact diagonalization up to iV = 36. 
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lx 



FIG. 9: Zigzag area on a lattice with T\ = l x u and T2 = l v v. 



TABLE II: Values of S[RK] - S[RK;p] on a zigzag area, 
divided by In 2. The sector p used in Fig. [l(] is indicated by 
*, and gives the best estimate in most cases. 







(5[RK] - 


S[RK;p])/ln2 




lx ( — ly) 


++ 


-+ 


+- 




4 


1.0024* 


0.8051 


1.4910 


0.8051 


6 


1.0315 


0.9248 


1.0315 


1.0212* 


8 


0.9944* 


1.0022 


1.0017 


1.0022 


10 


0.9981 


1.0028 


0.9981 


1.0011* 



Since the area is topologically non-trivial (it contains 
the incontractible cut Ai), the value of Sn depends on 
the choice of the ground state, even for large systems. 
We calculate the entanglement entropies on this area in 
the ground-states |RK) and |RK;p) on isotropic lattices 
l x = l y , and write them as 5[RK] and SplK;??] respec- 
tively. The results are plotted in Fig. [Io[ As anticipated, 
SV2 appears to be almost perfectly linear in l x (compared 
with the results of Fig. [|). Moreover, we observe that the 
topological constant 7 can be extracted in two different 
ways: a) by extrapolating (through a linear fit) S"[RK;p] 
at l % = 0" or b) by -7 ~ S[RK;p] - S[BK]. These two 
follow from the scaling forms 

S = ai l x for |RK), 
S = ail x ~j for |RK;p), 

where u\ is a non-universal constant. A_similar scaling 
was obtained rigorously by Hamma et al£3 for a "ladder" 
area in Kitaev's model on the square lattice. Here we 
confirmed that it holds accurately even in a system with 
a finite correlation length. The scaling forms ([ll|) provide 
an accurate way to calculate the topological constant 7 
even in relatively small systems. The condition (satisfied 
by QDM) is that topological sectors must be well defined 
and not mixed by the Hamiltonian, so that one can label 
the ground-states by their sectors. Computing 7 from the 
largest system (l x = l y = 10) gives our best estimate of 
the topological entanglement entropy S*[RK] — S[RK;p — 

] = 0.6939 (to be compared with In 2 = 0.6931); see 

Table [n]. 

As another application, we use the scaling forms jll| ) 
to evaluate the topological term —7 in the region v/t < 1. 




_2 1 1 1 1 1 1 1 

2 4 6 8 10 

lx 

FIG. 10: (color online) Entanglement entropies on a zigzag 
area at the RK point. The upper line is a linear fit to S[RK]. 
The lower one is a fit to 5*[RK;p = H — h] when l x = l y is 

multiple of 4 and S'fRKjjj = ] otherwise. The topological 

constant estimated from the latter fit is 7 = 0.73 ±0.01. This 
particular choice of p as a function of l x is motivated by the 
fact that, when v/t < 1, it corresponds to the ground-state 
sector. Explicit values of S^RK] — S[RK; p] are shown in Table 

We performed Lanczos diagonalization for lattices with 
l x = l y — 4 and l x — l y — 6. For v/t < 1, the 
ground-state lies in the sector p = ++ for l x = 4 and 

p = for l x = 6. We therefore compute the entropies 

S[p = ++; l x = 4] and S[p = ; l x = 6] on the zigzag 

areas and approximate the topological term —7 by a lin- 
ear extrapolation to u l x — 0" . In the thermodynamic 
limit, the constant term extracted in this way is expected 
to jump from — In 2 to a positive value, as in the case of 
Fig. g. However, the \/l2 x \/l2 VBC ordering is com- 
patible only with lattices where l x = l y is a multiple of 
6, and a linear relation Sq ~ a\l x + lnd holds only for 
such lattices.EZI The lattice with l x = l y = 4 is thus out of 
this scaling, and the present estimation of the constant 
term is invalid for the VBC phase. Still, it can be used 
in the liquid phase. The result is shown in Fig. |ll|. A 
value close to In 2 is recovered at the RK point but it de- 
creases smoothly when decreasing v/t. No clear signature 
of a transition out of the topological liquid can be seen. 
Larger system sizes are probably required to locate the 
transition with this method. The problem probably lies 
in a rapid increase of the dimer-dimer correlation length 
(and thus stronger finite-size effects) when moving away 
from the RK point in the direction of the VBC phase. 



IV. SUMMARY AND CONCLUSIONS 

The concept of topological entanglement entropy was 
recently introduced by KP and LW as a way to detect 
and characterize topological order from a ground-state 
wave-function. We have illustrated numerically how this 
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-0.2 
-0.3 
-0.4 



Zigzag area, ED (N=4x4, 6x6) 

extrapolation using two points: 
-T=(6S[++;l x =4]-4S[--;l x =6])/(6-4) 




FIG. 11: The topological term —7 as a function of v/t, esti- 
mated by using zigzag areas for l x — 4 and 6. 



approach works in the case of the Z2 liquid phase of the 
QDM on the triangular lattice. We found that, due to 
lattice discretization, the topological entropy 7 cannot be 
obtained from a direct fit to the scaling form S ~ qL- 7. 
Instead, it is necessary to combine the entropies on plu- 
ral areas to cancel out the boundary contributions, as 
suggested by KP and LW. In particular, for the KP con- 
struction, we clearly observed that in the large-area limit 
the topological entanglement entropy converges to — In 2 
expected for Z2 topological order. We also proposed a 
procedure to evaluate the topological entropy using a 
topologically non-trivial "zigzag" area, which gives an ac- 
curate value even in small systems. For a system of linear 
size l x = 10, the later method provided an estimate of 
the topological entanglement entropy 0.6939 at the RK 
point, in remarkable agreement with the expected value 
(In 2 = 0.6931). 

In addition to illustrating the concept of topological 
entanglement entropy in a "realistic" model, the present 
analysis offers an evidence of Z2 topological order in the 
QDM on the triangular lattice from a new perspective. 
Although the existence of topological degeneracy,li3 the 
analogy between this model and a Z 2 gauge theorye2l 
and the absence of any broken symmetryQ were already 
known, the present work confirms the Z2 structure in the 
ground-state wave-function itself. 
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APPENDIX A: REDUCED DENSITY MATRIX 
OF THE RK WAVE FUNCTION 

In this appendix, we derive a simple expression for the 
RDM of the RK wave function (ji|), and describe two 
methods for calculating it. A dimer configuration C on 
the entire system can be divided into the configurations 
on Q and Q: 



(Al) 



where £q (£q) denote the set of all the possible dimer 
configurations on fl (fi). Now we consider the inverse 
mapping: given c € £n and c £ £q, under what condition 
is (c, c) a physical configuration? This condition is given 
in terms of "occupied sites" of dimer configurations as 
follows. For a configuration c G £q, we define A(c) as the 
set of all the sites occupied by dimers in c, as shown in 
Fig. [l^. We similarly define A(c) for c £ £q. In order for 
(c, c) to be physical, a) A(c) and A(c) should not overlap 
with each other, and b) the sum of A(c) and A(c) should 
cover all the sites of the lattice. Then we can rewrite the 
wave function (Q) as 



IRK) = -= £ |c) 

A(c) U A(c) = X s 



(A2) 



where X s is the set of all the sites. If we list up all the 
possible A(c) and write them as A, (i — 1, 2, • • • ), we can 
divide the summation as 



|RK> 



^E E 14 E 

1 c e £q cefjj 

A(c) = A; A(c) = X s \ A; 



(A3) 



We introduce 

4 = {cef n |A(c) = A i } ) 

4 = {cefn|A(c) = X,\A i }, 

and we define normalized states on Q and Cl as 



1 



Ei c >> \4 



1 



Eic>- 



c£E' n 



Then we arrive at the Schmidt decomposition : 

141 • n 



|RK)=£v^l</4>r4>, with a 



(A4) 



(A5) 



(A6) 
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FIG. 12: Left (right): dimer configuration c (c) on Q (fi) 

and their "occupied sites" A(c) (A(c)), marked with circles 

(squares). In this picture, A(c) and A(c) are compatible and 

thus (c, c) is physical. 



The RDM of this state reads 

Pn = Tr|RK)(RK| = ^(^|RK)(RK|^) 

i 

= ^|^)A l (^|. 



(A7) 



This expression is already diagonal and the Aj's are the 
eigenvalues of pn. The entanglement entropy is then 
given by Sq = — J^AjlnAj. Since Aj's are expressed 
using the number of dimer coverings for a given set of oc- 
cupied sites, the task has been reduced to counting dimer 
coverings. This can be done by direct enumeration using 
a recursive algorithm or by a Pfaffian method. 

The Pfaffian method uses the fact that the number of 
dimer coverings is given by the the Pfaffian of an adja- 
cency matrix with appropriate signs (entries .are ±1 if 
the two sites are connected, and otherwise) .t3 Count- 
ing only configurations (c, c') such that A(c) = Aj and 
A(c) = X \ Aj can be done by removing some bonds of 
the lattice (setting to zero the corresponding matrix el- 
ement): if a site x belongs to Aj, there cannot be any 
dimer between x and a site y if (xy) is not a bond of 
Q. In the same way, any bond (xy) S involving a site 
x Aj must be switched off. The product |£q| • \£^\ is 
thus obtained Jxpm the Pfaffian of the modified adjacency 
matrix above. c3 It is clear that sites in the "bulk" of fl 
are necessarily included in all Aj (otherwise the number 
of configurations is zero) and those in the bulk of are 



necessarily excluded. There is a choice only for the sites 
in the vicinity of the boundary between Q and f2 (sites 
which arc both connected to bonds € fi and bonds ^ fi). 
The number of possible Aj therefore scales as 2 p ( n ) where 
P(f2) is the "perimeter" of ft. Since the calculation of 
each Pfaffian requires of the order of ~ iV 3 operations 
(see Ref. ^ for an explicit algorithm) the computer time 
required to obtain the RDM (and its spectrum) scales as 
~ A^ 3 • 2 p ( n ). This method is thus appropriate to study 
"small" areas in "large" systems. The results of Fig. 
for zigzag areas with l x = 8 and l x — 10 were obtained 
by this method. 

The direct enumeration algorithm searches and counts 
physical dimer configurations one by one for a given set 
of occupied sites. The enumeration is done separately for 

and £l, and the required time for each area is almost 
proportional to the number of dimer coverings, \£q\ or 
\£h\- Let AT(f2) be the number of sites in the "bulk" 

of n 



then \£a\ scales as 



a is 



a con- 



stant. Similarly, \£±\ ~ a N ( n \ Since we have 
possible Aj's, the total computation time adds up to ~ 
2 P(n)( a N(Q) + a N(Q)). With the extension of the area fi, 
the number of possible Aj increases, but counting dimer 
configurations get faster because ft shrinks. Thus this 
method is optimal for large areas in medium-size systems 
(here up to N — 64), being complementary to the Pfaf- 
fian method. One can reduce the time further by dividing 
fl or f2. Let us consider an annulus like in Fig. [| as f2, 
for example. Then can naturally be divided into inner 
(r < i?i n ) and outer (r > Rout) parts, denoted by u> and 
a/. Since f2 has two disconnected boundaries, with to and 
with a/, one can label A(c) by two numbers, i and j, cor- 
responding to the occupations around these boundaries. 
The dimer configurations on w and u>' can be counted 
separately for given i and j. The eigenvalues to calcu- 
late is therefore expressed as Ay = |££| • |£q | • 
Let P (P 1 ) be the "length" of the boundary between 
and lj (u)'). The required computation time becomes 
~ 2 P ■ a N ^ + 2 P+P ' ■ a N ^ + 2 P ' ■ a N ^'l By dividing 
areas, in general, one can reduce the time of counting 
configurations in this way, but the number of possible 
occupations at the boundaries increases. One needs to 
choose an efficient division, depending on the system and 
area sizes. 
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